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ON  K-LINE  AND  K  x  K  BLOCK  ITERATIVE  SCHEMES 
FOR  A  PROBLEM  ARISING  IN  3-D  ELLIPTIC  DIFFERENCE  EQUATIONS 

by 

Seymour  V.  Parter 
and 

Michael  Steuerwalt 
ABSTRACT 

Novel  computer  architectures  and  a  desire  to  solve  three- 
dimensional  problems  have  together  aroused  new  interest 
in  iterative  methods  for  computing  solutions  to  elliptic 
difference  equations.  Block  iterative  methods  are 
particularly  attractive  for  vector  machines ,  such  as  the 
CRAY-1.  Plane  iterative  schemes  reduce  a  three- 
dimensional  elliptic  system  to  two-dimensional  systems. 

We  analyze  the  convergence  rate  of  k-line  and  k  x  k 
block  iterative  methods  for  solving  these  two-dimensional 
problems. 

\ 

\ 

1.  INTRODUCTION 

Some  years  ago  there  was  great  interest  in  iterative  methods 
for  solving  elliptic  difference  equations:  see  [1,  12,  13,  21,  22, 
23,  24].  Recently  there  has  been  more  emphasis  on  finite  element 
equations  ([2,  4,  S,  20,  25])  and  direct  methods  of  solution 
C[7,  8,  11,  17,  18]).  Nevertheless,  in  practice,  particularly  in 
the  case  of  three-dimensional  problems  (see  [9,  11,  19]),  finite 
difference  equations  are  frequently  used  and  they  are  commonly 
solved  by  an  iterative  method,  usually  some  variant  of  the  SOR 
method.  Furthermore,  the  advent  of  new  computer  architectures, 
"vector  machines"  and  "parallel  processors,"  has  stimulated  a 
search  for  iterative  schemes  that  are  particularly  efficient  for 
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these  machines. 


•hub** 


In  this  report  we  examine  two  problems  that  arise  in  this 
way.  Consider  the  three-dimensional  model  problem 

(1.1)  -  Ah(3)  U  •  f  ,  (xit  Yy  tn)  6  ft 

Cl- 2)  u  -  g  ,  O4,  y^ ,  tn)  6  3ft 

where 

(1.3a)  (xi#  Yy  zn)  ■  (ih,  jh,  nh)  ,  0  <  i,j,n  <  P  ♦  1 

are  the  grid  points  in  the  unit  cube  Q  with 
(1.3b)  h  -  __ L-r  * 


and  Ah(3)  is  the  usual  7  point  discrete  approximation  to  the 
Laplace  operator,  given  by 


r 


d.4)  [Ah( 3)U]iJn 


Ui*l, j  ,n  *  2  Uijn  *  ui-l,j,n  + 
h2 

Ui.j*l,n  “  2  Uijn  *  Ui,j-l,n  + 

—  h2'-  — 

Ui,J,n+l  '  2  uijn  *  Ui,j,n-1 


Suppose  that  one  has  decided  to  use  a  "plane"  Iterative  method 
probably  block  SOR  where  each  block  is  the  set  of  unknowns 
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associated  with  a  plane  n  •  constant.  New  we  must  solve  the 
equations  in  each  plane.  With  V  -  UR,  these  equations  can  be 
written  as 


(1.6) 


6  Vij  ‘  Vi*lJ  '  Vi-l,j  ’  Vi,j*l  “  Vi,j-1  "  Fij» 

1  <  i.j  <P  . 


We  consider  two  block  iterative  schemes  for  the  solution  of  (1.6). 
These  are 

(i)  k-line  iterative  scheme:  In  this  scheme  each  block  con¬ 
sists  of  the  unknowns  associated  with  the  points  on  k 
consecutive  horizontal  lines.  These  blocks  are  indexed  by  a 
single  index,  say  s.  We  have 


(1.7) 


V  .  (y. 
s  lvi, 


k(s-l)+u 


♦u  ;  1  i  u  i  k)  * 


(ii)  k  x  k  block  iterative  scheme:  In  this  scheme  each  block 
consists  of  the  unknowns  Vjj  associated  with  the  points  in  a 
k  x  k  square.  It  is  easiest  to  describe  this  block  with  a 
double  index  (r,s).  The  (r,s)  block  is 

(i.8)  $rs  -  {vkCr.1)+i>k(s.1)+j  ;  1  1  i  kJ  • 

Because  each  of  these  block  structures  satisfies  block 
property  A  (see  (1,  22,  24])  the  spectral  radii  of  the  Gauss-Seidel 
method  and  the  SOR  method  (as  well  as  the  optimal  overrelaxation 
parameter  u^)  are  determined  by  the  spectral  radius  of  the  block 
Jacobi  scheme. 

From  the  analysis  of  the  corresponding  block  iterative  methods 
applied  to  the  two-dimensional  Poisson  problem  (see  [3,  14,  IS]), 
one  might  expect  that  these  spectral  radii  behave  like  1/k. 
Unfortunately,  this  is  not  the  case.  Zn  fact,  if  pCSk L)  and 
p (SkB)  are  the  spectTal  radii  of  the  block  Jacobi  iterative  methods 
based  on  the  blocks  (1.7)  and  Cl.®)  respectively,  applied  to  (1.6), 
we  find  that  both  p(SkL)  and  p(SkB)  have  nonzero  limits  as  k  *  •  ! 


In  Section  3  we  give  the  results  £or  the  case  k  •  1.  These 
results  are  obtained  immediately  by  tensor  product  methods. 

In  Section  4  we  begin  the  discussion  of  our  analysis  of  p(SkL) 
and  p(SkB).  In  particular  we  use  a  variant  of  "separation  of 
variables"  to  reduce  the  problem  to  a  one -dimensional  eigenvalue 
problem  involving  tridiagonal  matrices. 

In  Section  5  the  basic  estimates  are  obtained.  Loosely  speak¬ 
ing,  we  obtain:  if  P/k  >_  2  then 

(1=9)  p(SkL)  -  2  -  >/T ♦  0(h2)  «  .267949  ♦  0(h2)  . 

p(SkB)  -*•  1/2  (3  -  V5j  ♦  0 (h2)  -  .381966  ♦  0(h2)  . 

Two  sets  of  approximate  values  for  P(SkL)  and  P(SkB)  are  given  in 
the  following  tables.  The  results  are  correct  up  to  a  term  which 
is  0  (h2) . 


k 

Pi(SkL) 

P2 (SkL) 

1 

1/2  -  3ir2h2/8 

1/2  -  3ir2h 

Z/8 

2 

1/3  ♦  0(h) 

1/3  ♦  0(h) 

3 

.268541  ±  .185(-1) 

.286398  t 

.684(-3) 

4 

.267987  ±  ,483(-2) 

.272772  ± 

. 443(-4) 

6 

.267949  i  . 344(-3) 

.268293  t 

. 221 (-6) 

8 

.276949  ±  .247(-4) 

.267974  t 

.114(-8) 

12 

.267949  ±  . 127(-6) 

.267949  t 

.  302  (- 13) 

14 

.267949  t  . 913(*8) 

.267949 

18 

.267949  ±  . 471 (- 10) 

24 

.267949  t  . 1 7 8 ( - 1 3) 

27 

.267949 

Figure  1. 
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k 

Pi(SkB) 

P2(SkB) 

1 

2/3  -  r 

V/3 

2/3  -  *2 

h2/3 

2 

1/2  ♦  0(h) 

1/2  ♦  0(h) 

3 

.384848 

t  .51S(-1) 

.432468 

t  .390(-2) 

4 

.382296 

*  .187 (-1) 

.400477 

±  .477 (-2) 

6 

.381972 

t  . 266 (- 2) 

.384625 

±  .934(-5) 

8 

.381966 

±  •  387(- 3) 

.382353 

±  .196(-6) 

12  1 

.381966 

t  .824 (- S) 

.381974 

±  . 888 (- 10) 

14 

.381966 

±  .120(- 5) 

.381967 

±  .189 (-11) 

16 

.381966 

t  .175 (-6) 

.381966 

t  . 409 (- 13) 

18 

.381966 

±  .256 (-7) 

.381966 

24 

.381966 

±  . 795 (-10) 

30 

.381966 

t  .247(-12) 

36 

.381966 

Figure  2. 

tables 

the  shorthand  .185(-1) 

stands  for  .185  *  10 

Note  that  the  tolerances  for  P2  are  smaller  than  the  tolerances 
for  ,  even  though  the  estimates  Pj  appear  to  converge  more 
rapidly. 

As  we  indicated  above,  the  Jacobi  spectral  radius  determines 
the  spectral  radii  of  the  Gauss -Seidel  and  SOR  methods.  For  the 
Gauss-Seidel  scheme,  neglecting  0(h2)  terms, 

(1.10)  p(S»L-GS)  -  p (S-L) 2  -  .071797 

p(S-B-GS)  •  p(S-B)2  -  .145898  . 


The  optimal  w  for  the  SOR  method  is  given  by 


where  p  is  the  spectral  radius  of  the  Jacobi  method  (see  [22]), 
and  with  this  choice  of  w 

p(SOR)  -  <ub  -  1  . 


Hence 

(1.11)  p(S-L-SOR)  -  .018624 

p(S-B-SOR)  -  .039406  . 

Figures  1  and  2  show  that  in  a  plane  iterative  method  the  inner 
iterations  to  solve  (1.6)  need  not  use  a  very  large  k.  Indeed, 
at  k  *  8  the  spectral  radii  have  essentially  reached  their 
asymptotic  values. 

In  Section  6  we  describe  some  computational  results.  Finally, 
in  Section  7  we  comment  on  the  results  of  this  work. 

However,  one  important  comment  should  be  made  at  this  point. 

The  analysis  given  here  seems  to  be  very  special  and  possibly 
limited  to  model  problems.  On  the  other  hand,  standard  elliptic 
problems  can  now  be  handled  in  great  generality  by  the  methods  of 
£3,  13,  IS].  Unfortunately ,  it  is  not  difficult  to  see  that  those 
methods  do  not,  and  apparently  cannot,  apply  to  these  strongly 
diagonally  dominant  problems.  The  elliptic  problems  discussed  in 
the  earlier  works  are  regular  problems  while  the  system  (1.6) 
corresponds  a  two* dimensional  discrete  singular  perturbation  pro¬ 
blem 

(1.12)  -  h2  6h(2)  U  ♦  2  U  ■  F  . 

For  this  reason,  and  to  distinguish  from  the  notation  used  in 
(15],  we  designate  the  spectral  radii  as  p(SkL)  and  p(SkB). 

We  are  indebted  to  Dan  Boley,  Bill  Buzbee,  and  Molly  Mahaffy 
for  much  support  and  encouragement  during  the  evolution  of  this 
work.  Bill  Buzbee  aroused  our  Interest  in  the  problem  and  provided 
the  basic  support,  as  well  as  many  fruitful  discussions.  Dan 
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Boley  wrote  the  original  code  which  was  used  for  experimentation 
and  optimization  of  the  results  of  [3] .  Molly  Mahaffy  revised 
Can  Boley' s  code  and  carried  out  the  computations  described  here. 
These  computations  were  all  performed  on  the  CRAY-1  at  the  Los 
Alamos  Scientific  Laboratory. 

2.  Iterative  Methods 

In  this  section  we  describe  the  basic  iterative  schemes  of 
interest.  Consider  a  system  of  linear  equations 

(2.1)  AX  -  Y 

where  A  is  an  R  x  R  matrix.  Block  iterative  schemes  for  the 
solution  of  (2.1)  are  completely  described  by  describing  a  par¬ 
tition  of  the  R- vectors  into  blocks.  Specifically,  suppose  we 
imagine  all  R-vectors  X  partitioned  into  block  vectors  of  the  form 

(2.2)  X  -  (Xx,  X2 . Xr’)1 

where  each  X^  is  itself  an  Rj -vector  and  (of  course) 


(2.3) 


R  . 


Corresponding  to  this  partition  of  vectors  the  matrix  A  is  natu 
rally  partitioned  into  blocks 


where  each  Aj*  is  itself  an  R^  *  Rj  matrix.  In  particular,  each 
diagonal  block  A^  is  a  square  R^  *  R^  matrix.  The  block  Jacobi 
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iterative  scheme  associated  with  this  block  structure  is  now  given 
by 


(2.4) 


x.(v)  *  • 


The  corresponding  Gauss -Seidel  iterative  scheme  is 


C2.S) 


Ajj  xi 


rij  •  • £,  a3»  x*tv<1)  -  s  aj«  x»v)  *  y3  • 


j<s 


while  the  successive  overrelaxation  (SOR)  iterative  scheme  with 
overrelaxation  parameter  u  is 


(2.6)  A  xfv+15  -  -  u  £  A.  X^v+1)  -  w  L  Ajs  X 

jj  1  j>s  5  j<s  J 

♦  Ajj  (1  -  u)  X^v)  ♦  m  Yj  . 


(v) 

s 


In  every  case  we  have  a  splitting 


(2.7)  A  -  M  -  N  , 

and  the  iterative  scheme  takes  the  form 

(2.8)  M  xCv+1)  -  n  XCv)  ♦  Y  . 

In  particular,  for  the  block  Jacobi  scheme 

(2.9)  M  •  diagonal  (A^) . 


For  any  such  splitting,  let 

(2.10)  P  ■  Bax  {|X|  ;  det(XM  -  N)  ■  o}  . 


8 


It  is  well  known  ([6,  10,  22])  that  if  A  is  nonsingular  then  the 
iterates  xCv)  converge  to  the  unique  solution  X  of  (2.1)  inde¬ 
pendently  of  x(°)  if  and  only  if 


(2.11)  p  <  1  . 

Returning  to  our  problem  (1.6),  we  see  that  our  vector  X  is 
an  R-vector,  with  R  ■  P^,  corresponding  to  the  two-dimensional  grid 
vector  0n  *  V  .  Furthermore,  (1.7)  and  1.8)  define  two  distinct 
partitions  of  V. 

It  is  not  difficult  to  see  that  both  partitions  lead  to  Jacobi 
iterative  schemes  that  satisfy  block  property  A.  Hence  the  Gauss- 
Seidel  spectral  radius  and  the  SOR  spectral  radius,  as  well  as  the 
optimal  choice  of  w,  are  determined  by  the  spectral  radius  p  of 
the  Jacobi  iterative  scheme. 

The  problem  studied  in  this  report  is:  for  each  of  the  block 
Jacobi  schemes,  k-line  and  k  x  k  block,  determine  the  asymptotic 
behavior  of  p  as  P  •  (i.e.,  as  h  ♦  0). 

We  therefore  study  the  generalized  eigenvalue  problem 

(2.12)  X  M  V  *  N  V 

where  M  is  given  by  (2.9)  for  the  two  partitions  (1.7),  (1.8).  In 
addition  to  block  property  A  these  splittings  also  have  the  follow¬ 
ing  properties: 


(2.13a)  M  -  M*  -  M1 

and  M  is  positive  definite; 

(2.13b)  N  -  N*  •  N*; 

(2.14)  M  is  an  M-matrix  --  its  inverse  M  •  (ft^)  satisfies 
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Finally, 


N  »  (N^)  also  satisfies 
(2.15)  Nij  i  0  • 

Other  properties  of  M  and  N  will  be  developed  as  needed. 

3.  Estimating  p,  k  ■  1 

The  case  k  ■  1  is  easily  handled  by  the  method  of  tensor  pro¬ 
ducts:  see  [12,  23].  Here  we  merely  record  the  results. 

Let  p(SlL)  denote  the  spectral  radius  of  the  Jacobi  iterative 
scheme  based  on  a  block  partition  into  single  lines.  Then 

(3-D  p(SlL)  -  yi-fgg&rc  «  \  (1  -  I  h2)  • 

Let  p(Sl2B)  denote  the  spectral  radius  of  the  1*1  block 
Jacobi  iterative  scheme.  This  is  "point"  relaxation  and  we  have 

(3.2)  p(Sl2B)  -  4 . J*  »  |  (1  -  \  it2  h2)  . 

4.  Estimating  p,  k  >  2.  Preliminaries 

In  this  section  we  develop  some  general  properties  of  our 
particular  Jacobi  iterative  schemes.  We  use  these  properties  to 
reduce  the  problems  to  one-dimensional  eigenvalue  problems 
essentially  via  "separation  of  variables." 

Because  our  splittings  satisfy  block  property  A  we  know 
that  if  X  is  an  eigenvalue  of  (2.12),  so  is  -X  (see  [22,  23,  24]). 
M  is  positive  definite  and  N  is  symmetric,  so  all  the  eigenvalues 
of  (2.12)  are  real  (see  [6]),  and  thus  p  is  characterized  by 

<«•»  •  ■  w 
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4Li 


where  (  ,  )  denotes  the  usual  vector  inner  product 

(4.2)  (X,Y)  «  -  £  Xi  yt  . 

Let  k  >  2  be  a  fixed  integer  less  than  P.  Assume  that  P  is 
chosen  so  that  k  divides  P  --  that  is, 

(4.3)  P  »  k  Q 
where  Q  is  an  integer. 

A  unified  approach  to  our  two  problems  is  provided  by  the 
following  convenient  representation  of  N  in  each  case.  Consider 
the  one -dimensional  operator  fl  acting  on  vectors 

(4.4)  4  *  ($i »  $2 »  •••»  4p) 

as  follows:  for  1  <  s  <  Q-l  ,  0  <  j  <  k-1  , 

0  »  2  <  j  <  k-1 

^ks+1,  j  ■  0  , 

♦k5  ’  >  ’  1  • 

Let  Nx  and  be  the  two-dimensional  operators  which  act  on  grid 
vectors  V  ■  (V^j)  in  the  following  manner:  Nx  acts  on  V  only  in 
the  x  direction,  i.e.,  the  first  subscript,  and  in  that  direction 
Nx  acts  as  N.  Similarly  Ny  acts  on  V  only  in  the  y  direction. 

For  example,  for  1  <  s  <  Q*l,  O^i^k-1 

2  <  i  <  k-1  , 

i  -  0  , 

i  -  1  . 


(4.6) 


(Nx  V>ks*i,;J 


ks+1, j  * 
^ks,j  » 


(4.S)  (»  -  J 
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A  careful  examination  of  the  two  partitions  of  the  matrix  A 
yields  the  first  lemma. 

Lemma  4.1:  For  the  k-line  Jacobi  iterative  scheme 
(4.7a)  N  -  Ny  . 

For  the  k  *  k  block  Jacobi  iterative  scheme 
(4.7b)  N  “  Nx  +  Ny  * 

Lemma  4.2:  Let  p(SkL)  and  p(SkB)  be  the  spectral  radii  of  the 
k-line  and  k  *  k  block  Jacobi  iterative  schemes  respectively. 

Then 

(4.8a)  *  °(h)  -  p(SkL)  i  7  *  0(h)  » 

(4.8b)  *  0(h)  -  p(SkB)  -  T  +  0(h)  * 

Proof:  We  first  obtain  the  lower  bounds  of  (4.8).  Let  U  »  (U-) 
be  the  grid  vector 

(4.9)  •  sin  itih  •  sin  irjh 

Then 


However,  U  is  an  eigenvector  for  A  and 

(4.10)  AU  •  (4  -  2  cos  *h)U  ■  (2  ♦  0(h2))U  , 
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whence 


(4.10) 


(NU ,U) _  # 

(2  +  0(h2))  (U,U)  ♦  (NU ,U) 


A  direct  calculation  using  the  "smoothness"  of  U  and  the  form  of 
ff  (see  [3,  14,  IS])  shows  that 


(4.11) 


(Nx  U.U)  -  (2/k  ♦  0(h))  (U,U) 
(Ny  U.U)  -  (2/k  ♦  0(h))  (U,U) 


We  put  (4.7)  and  (4.11)  into  (4.10)  and  use  (4.1)  to  obtain  the 
lower  bounds  of  (4.8). 

Let  V  be  an  eigenvector  associated  with  the  eigenvalue  p.  Then 


p  M  V  -  N  V  . 


Subtracting  pNV  from  both  sides  gives 


(4.12a)  A  V  -  u  N  V 

where 

(4.12b)  u  ■  . 

Hence 

(4.13)  (AV,V)  •  u (NV, V) 


It  is  an  easy  matter  to  see,  using  the  explicit  eigenvalue  of  A 
or  the  Gerschgorin  theorem,  that 


(4.14a)  (AV,V)  >  2 (V ,V)  . 
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Moreover,  the  definition  of  N  shows  that 


(4.14b)  (Ny  V,V)  <  CV.V) 

(4.14c)  ([Nx  ♦  Ny]  V.v)  <  2  (V,V)  . 

Thus  (4.13)  and  (4.14a),  together  with  (4.14b)  or  (4.14c),  show 
that  2  <  ii  or  2  2 y,  respectively.  The  upper  bounds  of  (4.8) 

now  follow  at  once. 

Let 


(4.15) 


0  ’  b 


We  rewrite  (4.12a)  as 


(4.16) 


a  A  V  «  N  V 


Any  positive  eigenvalue  a  of  (4.16)  corresponds  to  a  positive 
eigenvalue  X  of  (2.12)  and  conversely  via  the  relationships 


(4.17a) 


•  r^T 


(4.17b) 


* 


Because  o  is  a  monotone  increasing  function  of  X  (and  conversely) 
we  seek  the  largest  positive  eigenvalue  a  of  (4.16).  But  A~*  is 
a  positive  matrix  and  N  f  C  is  a  nonnegative  matrix,  so  by  the 
Perron- Frobenius  theory  [22]  the  largest  positive  eigenvalue  of 

(4.16)  is  its  spectral  radius,  and  the  associated  eigenvector  V 
may  be  taken  nonnegative. 

Let  us  study  the  eigenvalue  problem  (4.16).  Because  A  is 
positive  definite  and  N  is  symmetric  there  is  a  complete  set  of 
eigenvectors  V^,  v  •  1,  2 . P2. 


We  attempt  to  apply  "separation  of  variables"  to  this  eigen¬ 
value  problem. 

Case  1:  The  k-line  scheme. 

For  each  n,  1  <  n  <  P,  let 

(4.18)  -  sin  irinh  ,  1  <  i,j  <  P  . 

Substitution  into  (4.16)  with  N  •  Ny  yields 

(4.19a)  a  An  4(n)  -  N  , 

where  An  is  the  tridiagonal  matrix 
(4.19b)  An  ■  [-1,  6-2  cos  irnh,  -1] 

Each  AR  is  positive  definite,  so  each  eigenvalue  problem  (4.19a) 
has  P  linearly  independent  eigenvectors,  say  4^(r),  r  •  1,  2, 

....  P.  Now  the  vectors  given  by 

(4.20)  »  sin  irinh  (r) 

2 

are  P  linearly  independent  eigenvectors  of  (4.16).  Hence  all 
the  eigenvalues  of  (4.16)  are  given  by  the  eigenvalues  of  the  P 
eigenvalue  problems  (4.19a),  n  •  1,  2,  ...,  P. 

We  therefore  seek  the  largest  positive  eigenvalue  o  of  the 
P  eigenvalue  problems  (4.19a).  However,  each  An  is  not  only 
positive  definite  but  also  an  M-matrix,  i.e..  A**  is  a  positive 
matrix.  Because  8  is  a  nonnegative  matrix,  the  largest  positive 
•igenvalue  of  (4.19a)  is  also  the  spectral  radius  of  that  problem. 
Moreover,  the  associated  eigenvector  4^  may  be  taken  nonnegative. 
Assume  that  has  been  done.  Then  both  V,  the  eigenvector  of  (4.16) 
associated  with  o,  and  4^,  the  eigenvector  of  (4.19a)  associated 
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with  o,  aust  be  nonnegative.  Therefore  the  representation  (4.20) 
shows  that  we  Bust  have 

n  ■  1 

We  summarize  these  facts  in  the  following  lemma. 

Lemma  4.3:  Consider  the  k-line  iterative  scheme  with  2  <  k  <  P 
and  (4.3)  holding.  Then 

p(SkL)  -  pS-5 

where  a  is  the  largest  positive  eigenvalue,  and  the  spectral  radius, 
of  the  eigenvalue  problem 

(4.21)  oA^  ♦  «  N  ♦ 

and  Ax  is  given  by  (4.19b)  with  n  ■  1. 

Case  2:  The  k  x  k  block  scheme. 

This  case  is  a  bit  more  complicated,  so  we  proceed  with  a 
more  formal  development  of  the  argument. 

Lemma  4.4:  Let  N  ■  Nx  ♦  Ny.  Consider  the  eigenvalue  problem 
(4.16)  Let  o  be  the  largest  positive  eigenvalue.  Then  o  is  a 
simple  eigenvalue:  there  is  only  one  linearly  independent 
eigenvector  associated  with  a.  Moreover,  the  associated  eigen¬ 
vector  may  be  taken  to  be  strictly  positive. 

Proof:  Consider  the  matrix  A'XN.  Because  A*1  is  a  positive 
matrix  and  N  is  a  nonnegative  matrix  not  identically  zero,  every 
column  of  A  *N  is  either  the  zero  vector  or  a  strictly  positive 
vector.  Let  T  be  the  permutation  matrix  which  collects  the 
positive  columns  into  the  first  r  columns,  so  that 

*11  0 

C4.22)  T*  A* XNT  -  ■  B  . 

*21  0 
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2 

Hare  is  in  r  *  r  positive  matrix  end  B2,  is  e  (P  -  r)  *  r 
positive  Mtrix.  Let 


be  an  eigenvector  of  B  with  associated  eigenvalue  X.  Then 


(4.23a) 

Bu  X  •  X  X 

(4.23b) 

B21  X  -  X  Y 

Thus  X  and  X  are  an  eigenvalue  and  associated  eigenvector  of  B^  . 
In  particular,  if  X  •  o  then,  because  B^  is  strictly  positive, 
there  is  only  one  linearly  independent  eigenvector  X  and  X  can  be 
taken  strictly  positive.  Since  (4.23b)  determines  Y  uniquely  in 
terms  of  X,  the  lemma  is  proven. 

We  are  now  ready  to  reduce  the  eigenvalue  problem  (4.16)  to 
a  one -dimensional  problem. 

Lemma  4.S:  Let  N  ■  Nx  +  Ny  consider  the  eigenvalue  problem 
(4.16).  Let  ?  be  the  largest  positive  eigenvalue.  Then  ^  is  also 
determined  as  the  largest  positive  eigenvalue  of  the  eigenvalue 
problem 

(4.24)  o  B  6  •  N  6 

where 

(4.25a)  ♦  •  $2 *  •••*  ^p) 

and  B  is  the  tridiagonal  matrix 
(4.25b)  B  -  (-1,  3,  -1] 
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Proof:  The  matrix  B  is  both  an  M-matrix  and  positive  definite. 
Therefore  the  eigenvalue  problem  (4.24)  has  P  linearly  independent 
eigenvectors.  Moreover,  if  oQ  is  the  largest  positive  eigenvalue 
the  associated  eigenvector  4°  may  be  taken  strictly  positive. 

Let  o  be  an  eigenvalue  of  (4.24)  and  let  4  be  an  associated 
eigenvector.  Let  the  gTid  vector  V  be  given  by 

(4.26)  Vi3  -  4*  ^  . 

Then 

o(AV)..  .  o  4i(B4)j  ♦  a  4j(B4)i  . 

We  apply  (4.24)  to  see  that 

9 (AV)  A j  -  4i(N4)j  ♦ 

-  (Ny  V)^  ♦  (Nx  V)^  -  (NV)^  . 

In  other  words,  the  grid  vector  V  is  an  eigenvector  of  (4.16)  with 
associated  eigenvalue  9.  In  particular,  if  0  »  oQ  and  4*4°  then 
the  grid  vector  V  is  not  only  an  eigenvector  of  (4.16),  it  is  a 
strictly  positive  eigenvector!  Hence,  by  virtue  of  lemma  4.4,  the 
V  so  obtained  is  an  eigenvector  of  (4.16)  associated  with  9,  the 
largest  positive  eigenvalue  of  (4,16).  This  proves  the  lemma. 

5.  Estimating  0,  k  >  3 

In  this  section  we  study  the  one -dimensional  eigenvalue  prob¬ 
lems  (4.21),  (4.24).  We  shall  first  reduce  the  problem  still 
further  by  eliminating  from  (4.21)  and  (4.24)  variables  corre¬ 
sponding  to  those  equations  in  which 

(N4)j  -  0 

In  order  to  do  this  we  require  a  specific  representation  of  the 
solution  of  tridiagonal  systems  of  equations. 
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TP— 


Let  k  >.  3  be  fixed. 

Lemma  S . 1 :  Consider  the  system  of  linear  equations 
(S.la)  -♦j.j  ♦  0^  -  •  0  ,  j  -  1,  2.  ....  k-2 

where  and  ^  are  given  and 


(5.1b) 

6  >  2  . 

Let  E^  ,  j  - 

0,  1#  ...» 

(5.2) 

E0  •  ° 

and  set 

(5.3) 

ak  -  Ej  E2 

bk  *  Ek-2 

Then 

(5.4a) 

♦l  -  ak  *0 

(5.4b) 

♦k-2  "  bk  1 

Furthermore 

C5.5) 

0  <  Ej  <  Ej 

,  k-2  be  generated  by  the  recursive  scheme 


5  m  1*  2,  . • .  |  k- 2 , 
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and  as  k  tends  to  • 


(5.6a)  bk  "  Ek-2  -  7  (*  '  ^  •  4  )  » 

(5. 6b)  a^  ♦  0  . 

Proof:  The  formulae  (5. 4a),  (5.4b)  are  obtained  from  the  well 
known  algorithm  for  solving  diagonally  dominant  tridiagonal  sys* 
terns;  see  (16,  10].  The  monotonicity  of  the  Ej  and  the  bound 
given  in  (5.5)  are  also  well  known  and  easily  established  by  induc¬ 
tion.  Finally,  if  k  ■*  •  then  E k_2  must  converge  to  E,,,,  a  solution 

of 


E2  -  6  E  ♦  1  ■  0  . 


Because  one  root  is  bigger  than  1  and  the  other  less  than  1,  the 
bound  (5.5)  implies  that  E.  must  be  the  smaller  root.  This  proves 
(5.6a).  The  proof  of  (,5.6b)  is  immediate  from  (5.5),  (5.6a)  and 
(5.3). 

Returning  to  the  eigenvalue  problems  (4.21)  and  (4.24),  let 
1  s  <  Q*2  and  consider  the  equations  satisfied  by  ♦jts+2*  ♦ks+3» 
...»  ♦ks+k-1*  *?or  tbe  COTresPonding  equations  we  have  (Sfo)j  •  0. 
Hence 


(5.7)  -  ♦ks+j-l  *  ^ks+j  *  *ks*j+l  "  0  * 

where  and  $k(s*l)  can  be  taken  as  known. 


i  -  2,  3, 


k-1 


In  these  equations 


(5.8a)  8*6-2  cos  irh 


for  the  eigenvalue  problem  (4.21)  and 
(5.8b)  6-3 
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for  the  eigenvalue  problem  (4.24).  in  either  case  we  use  lemma 

5.1  together  with  the  equations  numbered  ks  for  s  •  2,  3 . Q-l 

and  the  equations  numbered  ks  ♦  1  for  s  ■  1,  2,  . ..,  Q-2  to 

eliminate  ,  j  •  2,  5,  . . . ,  k-1  and  s  -  1,  2 . Q-2.  For 

example,  the  ks  equation  is 

0  ["  ♦ks-l  +  8*ks  ’  ♦ks+l]  *  *ks  +  l 

If  2  <  s  £  Q-l,  then  with  the  appropriate  choice  of  a^  and  bk  we 
have 

♦ks-l  *  bk  0ks  *  ak  ♦k(s-l)+l  • 

We  thus  obtain  foT  2  <_  s  £  Q-l 

(5‘9)  0  [’  ak  ♦k(s-l)+l  +  (S  ‘  V  ♦ks  *  ♦ks+l]  "  ♦ks+l  * 

Similarly,  the  (ks  +  1)  equation  is 

0  ["  ♦ks  +  6  ♦ks+l  ‘  ♦ks+2]  "  ♦ks  * 

If  1  <  s  <  Q-2,  then  with  the  appropriate  choice  of  ak  and  b^ 
we  have 


♦ks+2  *  ak  ♦k(s+l)  +  bk  ♦ks+l  * 
Thus  we  have,  for  1  £  s  <  Q-2, 


(5.10) 


o  [-  ♦ks  ♦  (S  -  t>k^ks*l  *  ak  ♦k(s+l)]  "  4ks  • 
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Now  we  must  eliminate  *2,  ....  #k-1  and  ♦k(Q-i}+2. 

*k(Q-l)+3»  •••»  *kQ*  In  these  cases  we  have  a  system  of  k  •  1 
unknowns.  However,  the  procedure  is  exactly  the  same.  We  collect 
our  results  in 

Lemma  5.2.  Let  o  f  0  and  4  be  an  eigenvalue  and  eigenvector 
respectively  of  (4.21)  or  (4.24).  Let  ak,  bk,  and  bk+1  be  given 
by  (5.3)  with  the  appropriate  choice  of  6,  either  (5.8a)  or  (5.8b). 
Let 

(5.11a)  52s-1  ‘  *ks  »  s  •  1,  2,  ...»  Q-l  , 

CS.Hb)  ?2s  *  ^ks+1  »  s  -  1,  2 . Q-l  . 

Let 

(5.11c)  y  •  l/o  ,  y  ■  1  ♦  u 

Then  y  and  £lt  ....  ^(Q-l)  satisfy  the  homogeneous  system 
(5.12a)  (0  -  bk+1)51  •  Y  52  ■  0 

(5.12b)  *  y  ♦  (8-  bk)C2  -  ak  ■  0 

(5.12c)  -  ak  52s  *  (8  -  bfc)€2s+]>  *  TC 2s+2  “  °* 

s  ■  1,  2,  . . . ,  Q-3 

(5 . 12d)  -  y  52$+1  ♦  (8  -  bk)52s4>2  -  akC2s+3  -  0, 

s  *  1 i  2,  • • • j  Q-3 


(S.12*) 


“  *k  *2(Q-2)  *  *  bk^2(Q-2)+l  ‘  y  *2(Q-1)  *  0 


(5.12f) 


Y  52(Q-2)+l  *  “  bk+l)52CQ-l)  "  0 


Moreover,  let  Uq  be  the  smallest  positive  number  for  which  the 
system  (S.12)  possesses  a  nontrivial  solution;  then  Yg  >  1  and 


(5.13) 


P 


i  .  i_  . 

t0 


Proof:  It  is  only  necessary  to  verify  the  characterization  of 
or  Yq  given  by  (5.13).  However,  this  is  an  immediate  consequence 
of  our  earlier  characterization  of  a  as  the  largest  eigenvalue  of 
(4.21)  or  (4.24). 

Corollary :  Consider  equations  (5.12)  with  bj{+1  replaced  by  bk< 
Let  y  be  the  smallest  positive  value  >  1  of  y  so  that  these  modi¬ 
fied  equations  have  a  nontrivial  solution.  Also,  consider  the 
system  (5.12)  with  b^  replaced  by  bj^j.  Let  Y  be  the  smallest 
positive  value  of  y  so  that  these  modified  equations  have  a  non¬ 
trivial  solution.  Then 


(5.14)  y  <  y0  <  y  . 

Proof :  Because 

bk  i  bk.l 


we  see  that 


When  y  is  a  small  positive  number  the  tridiagonal  matrix  of  (5.12) 
is  positive  definite.  If  we  increase  some  diagonal  elements  (for 
instance,  replace  g  -  bj^  by  B  -  bk)  we  increase  all  the  eigen- 
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values.  Thus  Increasing  some  diagonal  elements  means  we  must  raise 
Y  to  make  the  system  singular.  Therefore 

y0  <  y  • 

Similarly,  replacing  by  bj^  makes  the  diagonal  smaller  and  y 
need  not  be  so  large.  Thus  we  have  proven  (5.14). 

Having  es-tablished  this  corollary,  we  proceed  to  estimate  the 
quantities  y,  y  . 

We  first  rearrange  the  matrices  of  interest. 

Lemma  5.3:  Let  A  be  the  symmetric  tridiagonal  matrix  of  order 
2 (Q  -  1)  of  the  form 


That  is,  the  diagonal  of  A  is  a  constant,  g.  The  superdiagonal  is 
alternately  -y,  -a^,  etc. 

Let  T  be  the  permutation  matrix  corresponding  to  the  permuta¬ 
tion 

(5,15)  (2j  -  1)  -  i  ,  2)  -  (Q  -  1)  ♦  )  , 

so  that,  letting  e^  denote  the  jth  unit  vector, 
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at 


be  a  nontrivial  null  vector  of  5,  so  that 
C5.19)  5  V  -  0  . 

Then  C5 . 19 )  is  equivalent  to  the  pair  of  equations 

(S. 20a)  E*  E  Y  -  I2  Y  , 

(5.20b)  E  E*  X  •  8 2  X  . 

A  computation  shows  that 


»'*.  "  ►  1 


t  2  2 

That  is,  EE  is  the  constant  tridiagonal  matrix  [a^y,  a^  ♦  y  , 
except  that  the  ((Q  -  1) ,  (Q  -  1))  term  is  y2  rather  than  ak  + 
Lemma  5.4:  Let  y  and  y  be  defined  as  in  the  corollary  to  lemma 
5.2.  Then 


(5.22) 


6  ■  bk.i  • 

8  •  bk  -  ak(l  •  0(h2))  <  8  -  bk  * 


and  from  (5.13)  we  have 


(5.23)  1  <  _ 1 _  <p 

6  *  bk  +  ak  6  -  bk  -  ak  ♦  ak  0(h2) 

<  1 _  . 

6  '  bk+l  '  ak 

Proof:  Consider  the  case  of  y.  Then  8*8-  b^+1.  From  (5.20a) 
we  see  that  if  y  *  y  then  I2  is  an  eigenvalue  of  EtE.  Thus  by 
the  Gerschgorin  Theorem  [6] 

62  <  aj|  ♦  y2  ♦  2ak  y  •  (a^  +  y)2  . 

This  establishes  the  left  hand  inequality  of  (5.22). 

Consider  the  case  of  y.  Then  8  ■  8  -  b^.  For  y  <  y  the  matrix 
5  is  positive  definite.  Therefore  at  y  »  y  the  smallest  eigen¬ 
value  of  B  is  zero.  Thus  -I  is  the  smallest  eigenvalue  of 

0  E 
E*  0  . 

m 

If  ti  is  in  eigenvalue  of  Bfl,  so  is  -q.  Moreover,  n  is  an  eigen- 

*  W  2  t 

value  of  Bq  if  and  only  if  ^  is  an  eigenvalue  of  E  E.  We 
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conclude  when  y  •  Y  that  I2  is  the  largest  eigenvalue  of  E*E. 

Let  ^  ^  >  •••  >  Hq.j  be  the  eigenvalues  of  E*E.  Fro*  the 

inclusion  theorem  [6,  p.  149}  and  the  known  eigenvalues  of 

[  1,  0,  1  } 

we  see  that 

v  cos  r^r  -  Vi  • 1  i  J  i  o-2 


Thus 


This  proves  the  right  hand  inequality  of  (5.22). 

Remark:  It  is  important  to  observe  that  the  coefficient  of  the 
0(h2)  term  includes  afc.  Since  ak  0  rapidly,  this  term  is  truly 
negligible. 

Corollary:  Let  k  *►  •  .  Then 
(5.24)  p  -  \  |*  -  h2  -  4 

Proof:  Using  (5.6a)  and  (5.6b)  we  see  that 
P  *  |  \  b  is  k  ♦  •  . 
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However,  from  the  basic  equation  (5.2)  we  see  that 


0 


1 

7 


■  yp~r 


Theorem  5.1;  Consider  the  k-line  Jacobi  iterative  scheme  where 
k  divides  P  and  k  <  P.  Let  p(SkL)  denote  the  spectral  radius  of 
this  scheme.  Then  the  results  shown  in  Figure  1  are  correct  up 
to  a  term  which  is  0(h2) . 

Proof:  The  result  for  k  ■  1  follows  from  (3.1).  The  result  for 
k  *  2  follows  from  (4.8a).  The  results  for  k  ^  3  were  obtained 
from  a  computation  based  on  (5.23)  with  8*4.  The  column 
was  computed  with  the  coarse  lower  bound  of  (5.23),  and  p2  with 
the  fine  lower  bound. 

Theorem  5.2:  Consider  the  k  *  k  block  Jacobi  iterative  scheme 
where  k  divides  P  and  k  <  P.  Let  p(SkB)  denote  the  spectral  radius 
of  this  scheme.  Then  the  results  shown  in  Figure  2  are  correct. 
Proof:  The  result  for  k  -  1  follows  from  (3.2).  The  result  for 
k  •  2  follows  from  (4.8b).  The  results  for  k  >  3  were  obtained 
from  a  computation  based  on  (5.23)  with  0*3. 

6.  Computational  Results 

Using  codes  originally  prepared  by  D.  L.  Boley  (see  [3]),  Molly 
Mahaffy  computed  approximate  spectral  radii  for  the  Gauss -Seidel 
iterative  scheme  using  k  lines  and  k  x  k  blocks.  These  spectral 
radii  were  computed  by  the  power  method.  The  Gauss -Seidel  method 
was  chosen  because  the  general  theory  shows  that  the  Jacobi 
iterative  scheme  has  both  p  and  -p  as  eigenvalues  while  the  Gauss- 
Seidel  scheme  has  a  simple  eigenvalue  on  the  spectral  radius. 
Furthermore,  we  also  have 


(6.1) 


p (Gauss- Seidel) 


[p (Jacobi)]  2 


The  results  are  contained  in  the  following  tables.  In  all  cases 
P  »  128.  As  with  Figures  1  and  2,  the  columns  p2  and  p2  were 
computed  using  the  coarser  and  finer  bounds  of  (S. 23) ,  respectively. 


k-line  Iterative  Scheme 


k 

4 

8 

16 

32 

64 


p2  (Theory) 

P2  (Theory) 

p2 (Computed) 

.071840  t  .259(-2) 

.074404  ±  ,242(-4) 

.07662 

.071797  t  .132  (-4) 

.071810  t  .608 (-9) 

.07167 

.071797  t  . 351(-9) 

.071797 

.07164 

.071797 

.071797 

.07164 

.071797 

.071797 

.07164 

Figure  3. 


k  x  k  Block  Iterative  Scheme 


k  p2  (Theory) 

4  .146498  1  .143  (-1) 

8  . 145898  ±  .296(-3) 

16  .145898  ±  .134  (-6) 

32  .145898  t  .271(-13) 

64  .145898 


p2  (Theory)  p2  (Computed) 


160382 

t  .382 (-3) 

.16620 

146194 

±  .150 (-6) 

.14631 

145898 

t  .311  (-13) 

.14S90 

145898 

.14590 

145898 

.14S90 

Figure  4. 


7.  Comments 

It  is  of  interest  to  observe  that  the  results  of  Theorem  5.1 
and  Theorem  5.2  do  not  require  that  k/P  *  0.  In  fact,  those  results 
are  valid  as  long  as  k  <  P.  This  is  clearly  demonstrated  in  the 
computational  results.  For  example,  in  Figs.  3  and  4  we  see  that 
with  P  ■  128  and  k  ■  64,  the  largest  acceptable  value,  the  value 
of  p2  is  essentially  given  by  the  asymptotic  value  (k  ♦  •).  These 
results  seem  to  contradict  the  intuitive  feeling  that,  were  M  to 
include  a  much  larger  part  of  A,  the  spectral  radius  would  be  much 
smaller. 
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